#!/usr/bin/env perl

# LKUtils.pm: utility subroutine library
# Author: Lee Katz <Lkatz@cdc.gov>
# Some methods are taken from Andrey Kislyuk's AKUtils

package LKUtils;
require 5.12.0;
my ($VERSION) = ('$Id$' =~ /,v\s+(\d+\S+)/o);

use strict;
use warnings;
use List::Util qw(sum reduce shuffle);
#use List::Util qw(min max sum reduce shuffle);
use File::Basename qw/fileparse/;
use File::Spec;
use File::Copy;
use File::Temp ('tempdir');
use Data::Dumper; # useful for printing deep structures, e.g. to print %hash: print Data::Dumper->new([\%hash],[qw(hash)])->Indent(3)->Quotekeys(0)->Dump;
use Bio::Perl;
use URI::Escape;
use Time::Local;
use Config::General;
use threads;
use Thread::Queue;

use Exporter;
our @ISA = "Exporter";
our @methods = qw(logmsg mktempdir loadData loadCSV argmax filesDiffer fullPathToExec md5sum isFasta getNumCPUs loadConfig trim Shuffle array_unique seqHash file_put_contents file_get_contents getGff getGffHash readCrunch script_duration);
our %EXPORT_TAGS = (all => [@methods]);
Exporter::export_ok_tags('all');

local $SIG{'__DIE__'} = sub { my $e = $_[0]; $e =~ s/(at [^\s]+? line \d+\.$)/\nStopped $1/; die("$0: ".(caller(1))[3].": ".$e); };

END {

}

our $startTime=time();
our $fastaOut=Bio::SeqIO->new(-format=>"fasta");

sub script_duration{
  my $seconds=time()-$startTime;
  return $seconds;
}

sub array_unique{
  my %h;
  map { $h{$_}++ == 0 ? $_ : () } @_;
}

sub Shuffle{
  my $inputarray = shift;
  my $i;
  for ($i = @$inputarray; --$i;){
    my $j = int rand($i+1);
    next if $i==$j;
    my $tmp=$$inputarray[$i];
    $$inputarray[$i]=$$inputarray[$j];
    $$inputarray[$j]=$tmp;
  }
  return @$inputarray;
}

sub trim{
  my $str=shift;
  $str=~s/^\s+|\s+$//g;
  return $str;
}

sub phredToErr{
  my($p,$settings)=@_;
  if(ref($p) ne 'ARRAY'){
    $p=[split(//,$p)];
  }
  chomp(@$p);
  my @p=();
  for(@$p){
    next if(/^\s*$/);
    my $acc=exp( (ord($_) - 33) / ( -10 * log(10) ) );
    my $err=1-$acc;
    die "Error: there was a phred score $_ that results in a negative error rate $err\n" if($err<0);
    push(@p,$err);
  }
  return @p;
}

sub errToPhred{
  my($p,$settings)=@_;
  my @p=@$p;
  for(@$p){
    push(@p,
      -10*log($_)
    );
  }
  return @p;
}

###################################################
# functions below this line are taken or modified from AKUtils#
###################################################
sub logmsg{
  my ($msg,$settings)=@_;
  $|++;
  my $file=(caller(1))[1] || "unknownFile";
  $file=fileparse($file)."::";
  my $function=(caller(1))[3] || "main";
  my($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst)=localtime(time);
  my $time=timelocal($sec,$min,$hour,$mday,$mon,$year);
  my $ymd=join("-",$year,$mon,$mday);
  my $hms=join("-",$hour,$min,$sec);
  my $datetime=join(" ",$ymd,$hms);
  print "$datetime $file$function: $msg\n";
  $|--;
}


sub mktempdir(;$) {
	my ($settings) = @_;
	my $tempdir_path = File::Spec->join(File::Spec->tmpdir(), (split("::",(caller(1))[3]))[1].".$$.XXXXX");
	my $tempdir = tempdir($tempdir_path, CLEANUP => !($$settings{keep}));
	return $tempdir;
}

# Loads data from file using Data::Dumper
sub loadData($) {
	my ($file) = @_;
	open(IN, '<', $file) or die("Unable to open file $file for reading: $!");
	my $ref;
	{	no strict 'vars';
		local $/;
		my $d = <IN>;
		$ref = eval $d;
		use strict 'vars';
	}
	close IN;
	return $ref;
}

# Assumes the keys are unique. TODO: keyless load
# Expects header line, keys on column 0 unless otherwise specified (by column index or name)
# TODO: sanity checks
sub loadCSV($;$) {
	my ($file, $key) = @_;
	my %table; my $key_col;
	open(FH, '<', $file) or die("Could not open file $file for reading: ".$!);
	my $header = <FH>; chomp $header;
	my @fieldnames = split(/,/, $header);
	my %f; for (@fieldnames) { die("Empty or duplicate header field") if $f{$_} or $_ eq ''; $f{$_} = 1; }
	$key = 0 if not defined $key;
	for (0..$#fieldnames) { $key_col = $_ if $key eq $fieldnames[$_]; }
	$key_col = $key if $key =~ /^\d+$/;
	while (<FH>) {
		chomp;
		my @fields = split /,/;
		for (0..$#fieldnames) {
			$table{$fields[$key_col]}->{$fieldnames[$_]} = $fields[$_];
		}
	}
	close FH;
	return(\%table, \@fieldnames);
}

sub argmax {
	my ($hash) = @_;

	#my $bestval = -1e999;
  my($best,$bestval)=each(%$hash); # LK, untested
	while (my($key,$val) = each(%$hash)) {
		if ($val > $bestval) {
			$bestval = $val;
			$best = $key;
		}
	}
	return $best;
}

sub filesDiffer($$) {
# TODO: use Algorithm::Diff
# @f1=<FH>; $diff = Algorithm::Diff->new(\@f1, \@f2); next   if  $diff->Same();
	my ($file1, $file2) = @_;
	open(FH1, '<', $file1) or die("Could not open file $file1 for reading: ".$!);
	open(FH2, '<', $file2) or die("Could not open file $file2 for reading: ".$!);
	while (1) {
		my $line1 = <FH1>;
		my $line2 = <FH2>;
		last unless defined $line1 and defined $line2;
		return 1 if $line1 ne $line2;
	}
	return 1 if <FH1> or <FH2>;
	return 0;
}

# If argument is an executable in the current path, returns the full path to it, otherwise returns undef
sub fullPathToExec($) {
	my ($executable) = @_;
	my $fullpath;
	for ("", split(/:/, $ENV{PATH})) {
    my $possiblePath=$_."/".$executable;
		if (-x $possiblePath && ! -d $possiblePath) { 
      $fullpath = File::Spec->rel2abs($possiblePath); 
      last; 
    }
	}
	warn("Error finding full path to file ($executable)") unless (-x $fullpath); 
	return $fullpath;
}

sub md5sum($) {
	my ($filename) = @_;
	my $md5sum;
	die("$filename: file not found") unless -f $filename;
	if (`which md5sum 2>/dev/null`) {
		$md5sum = `md5sum $filename`; $md5sum =~ /^([^\s]+)/; $md5sum = $1;
	} else {
		require Digest;
		open(FH, '<', $filename) or die("Could not open file $filename for reading: ".$!);
		my $md5obj = Digest->new('MD5');
		$md5obj->addfile(\*FH);
		$md5sum = $md5obj->hexdigest();
		close FH;
	}
	return $md5sum;
}

# TODO: this needs more work
sub isFasta($) {
	my ($file) = @_;
	my $char;
	open(FH, "< $file") or die("Could not open file $file for reading: ".$!);
	read(FH, $char, 1); close FH;
	if ($char eq ">") { return 1; } else { return 0; }
}

sub getNumCPUs() {
	my $num_cpus;
	open(IN, '<', '/proc/cpuinfo'); while (<IN>) { /processor\s*\:\s*\d+/ or next; $num_cpus++; } close IN;
	return $num_cpus || 1;
}

# Load configuration values from config files.
# Usage: $settings = loadConfig($default_settings);
# By default, load 3 files in succession: bindir/appnamerc, /etc/appnamerc, ~/.appnamerc (if they exist)
# Override this with: $settings = loadConfig({config_files => [file1, file2]});
sub loadConfigAK($) {
	my ($settings) = @_;

	$0 = fileparse($0) if $0 =~ /\//;
	$$settings{appname} ||= $0;
	$$settings{app_config_dir} ||= "$FindBin::RealBin/../conf";
	$$settings{system_config_dir} ||= "/etc";
	$$settings{user_config_dir} ||= $ENV{HOME};
	$$settings{config_files} ||= ["$FindBin::RealBin/$$settings{appname}rc",
		"$$settings{app_config_dir}/$$settings{appname}rc",
        "$$settings{system_config_dir}/$$settings{appname}rc",
		"$$settings{user_config_dir}/.$$settings{appname}rc"];

	my @config_files = @{$$settings{config_files}}; # avoid ridiculousness from config_files being set in a config file

	foreach my $file (@config_files) {
		open(IN, '<', $file) or next;
		my $i;
		while (<IN>) {
			$i++;
			chomp;
			next if /^\s*\#/ or /^$/;
			/^\s*([^\=]+?)\s*=\s*([^\=]+?)\s*$/
				or die("Invalid configuration syntax on line $i of file $file. Expected variable=value syntax");
			my ($key, $value) = ($1, $2);
			if ($value =~ /^"/) { $value =~ s/^"//; $value =~ s/"$//; }
			$$settings{$key} = $value;
		}
		close IN;
	}

	return $settings;
}

# makes a hash of sequences, given a fasta file
# settings
#  seqHash_extra_keys => try to make redundant entries as a catch-all
sub seqHash{
  my ($filename,$settings)=@_;
  my($seqio,%seq);
  if(-e $filename){
    $seqio=Bio::SeqIO->new(-file=>$filename,-format=>"fasta");
  }
  # else it's not a file.  It could be a sequence object.
  elsif(ref($filename)=~/seq/i) {
    die ref($filename)." is not supported yet for seqhash";
  }
  else{
    die "Cannot determine the type of sequence, type ".ref($filename)."\n($filename)\n";
  }
  while(my $seq=$seqio->next_seq){
    my $id=$seq->id;
    $seq{$id}=$seq;

    if($$settings{seqHash_extra_keys}){
      # make a second key without genome info and zero padding
      my $id2=$id;
      $id2=~s/^.+?_0*//;
      $seq{$id2}=$seq;

      $id2=$id;
      $id2=~s/^(lcl|gi)\|//;
      $seq{$id2}=$seq;
    }
  }
  return %seq;
}

# reads a blast file
sub readBlast{
  my($blast,$settings)=@_;
  my @return;
  open(FILE,$blast) or die "Could not open blast file $blast because $!";
  while(<FILE>){
    chomp;
    my %t=(line=>$_);
    ($t{query},$t{subject},$t{identity},$t{length},$t{mismatches},$t{gaps},$t{qStart},$t{qStop},$t{sStart},$t{sStop},$t{evalue},$t{score})=split(/\s+/);
    
    # is it a +1/-1 match?
    if($t{qStart} < $t{qStop} && $t{sStart} < $t{sStop}){
      $t{strand}=1;
    }
    elsif($t{qStart} > $t{qStop} && $t{sStart} > $t{sStop}){
      $t{strand}=1;
    }
    elsif($t{qStart} > $t{qStop} && $t{sStart} < $t{sStop}){
      $t{strand}=-1;
    }
    elsif($t{qStart} < $t{qStop} && $t{sStart} > $t{sStop}){
      $t{strand}=-1;
    }

    push(@return,\%t);
  }
  close FILE;
  return @return;
}

# reads a crunch file
sub readCrunch{
  my($crunch,$settings)=@_;
  my @return;
  open(FILE,$crunch) or die "Could not open crunch file $crunch because $!";
  while(<FILE>){
    my %t;
    ($t{score},$t{identity},$t{qStart},$t{qEnd},$t{query},$t{sStart},$t{sEnd},$t{subject})=split(/\s+/);
    push(@return,\%t);
  }
  close FILE;
  return @return;
}

sub fastaOut{
  return $fastaOut;
}

sub file_put_contents{
  my($filename,$contents)=@_;
  $contents=join("",@$contents) if(ref($contents) eq 'ARRAY');
  open(FILE,">$filename") or die("Could not open file $filename for writing because $!\n");
  print FILE $contents;
  close FILE;
  return (-s $filename); # return bytes of the file
}
sub file_get_contents{
  my($filename)=@_;
  local $/; undef $/;
  open(FILE,$filename) or die("Could not open file $filename for reading because $!\n");
  my $str=<FILE>;
  close FILE;
  if(wantarray){
    return split(/\r*\n/,$str);
  }
  return $str;
}
# settings can have keys
#   splitSubattributes (default 1)
sub getGff{
  my ($file,$settings)=shift;
  my $gff;
  open(GFF,"<",$file) or die("Could not open the gff file $file because $!\n");
  my $numEntries=0;
  while(<GFF>){
    next if(/^#/);
    chomp;
    my $gffLine=LKUtils::readGffLine($_,$settings);
    push(@$gff,$gffLine);
    $numEntries++;
  }
  close GFF;
  
  return $gff;
}
# assumes unique contig/start key
# settings can have keys
#   splitSubattributes (default 1)
sub getGffHash{
  my ($file,$settings)=shift;
  my $gff={};
  #`wc -l $file`; # somehow helps stall until the full file is finished being written
  open(GFF,"<",$file) or die("Could not open the gff file $file because $!\n");
  my $numEntries=0;
  while(<GFF>){
    next if(/^#/);
    chomp;
    my $gffLine=LKUtils::readGffLine($_,$settings);
    $$gff{$$gffLine{seqname}}{$$gffLine{start}}=$gffLine;
    my $contigNumber=$$gffLine{seqname};
    $contigNumber=~s/contig0*//i;
    $$gff{$contigNumber}{$$gffLine{start}}=$gffLine;
    $numEntries++;
  }
  close GFF;
  
  return $gff;
}


# called by getGff, getGffHash
sub readGffLine{
  my($line,$settings)=@_;
  my $splitsubAttributes=$$settings{splitSubattributes}||1;
  my @t=split(/\t/,$line);
  my %gffLine=('seqname'=>$t[0],'source'=>$t[1],'feature'=>$t[2],'start'=>$t[3],'stop'=>$t[4],'score'=>$t[5],'strand'=>$t[6],'frame'=>$t[7]);
  #parse the attributes
  my $attribute=$t[8];
  my @attribute=split(/;/,$attribute);
  foreach(@attribute){
    my($key,$value)=split /=/;
    $key=uri_unescape($key);
    $value=uri_unescape($value) || $key;
    $gffLine{attribute}{$key}=$value;
    if($splitsubAttributes){
      my @tmp=split(/,/,$value);
      for my $t (@tmp){
        my($k,$v)=split(/:/,$t);
        next if(!defined($v));
        $gffLine{attribute}{$key."Hash"}{$k}=$v;
      }
    }
  }
  return \%gffLine;
}

# Load configuration values from config files.
# Usage: $settings = loadConfig($default_settings);
# By default, load 4 files in succession: bindir/appnamerc, /etc/appnamerc, ~/.appnamerc (if they exist), ./appnamerc
# Override this with: $settings = loadConfig({config_files => [file1, file2]});
sub loadConfig($) {
  my ($settings) = @_;
  my $realBin=$FindBin::RealBin || '.';

  $0 = fileparse($0) if $0 =~ /\//;
  $$settings{appname} ||= $0;

  $$settings{app_config_dir} ||= "$realBin/../conf";
  $$settings{system_config_dir} ||= "/etc";
  $$settings{user_config_dir} ||= $ENV{HOME};
  $$settings{environment_config_file} ||= "./$$settings{appname}rc" || "./.$$settings{appname}rc" || "";
  $$settings{config_files} ||= ["$realBin/$$settings{appname}rc",
    "$$settings{app_config_dir}/$$settings{appname}rc",
        "$$settings{system_config_dir}/$$settings{appname}rc",
    "$$settings{user_config_dir}/.$$settings{appname}rc",
    $$settings{environment_config_file}];

  # some standard settings to load before the config files are loaded
  $$settings{numcpus}||=LKUtils::getNumCPUs();

  my @config_files = @{$$settings{config_files}}; # avoid ridiculousness from config_files being set in a config file

  foreach my $file (@config_files) {
    next if(!-e $file);
    my %settings=Config::General::ParseConfig($file);
    %$settings=(%$settings,%settings);
  }

  return $settings;
}

#################
### analysis
#################
# TODO if wantarray, then provide a null Fst
sub Fst{
  my($seq1,$seq2,$settings)=@_;

  my %settings;
  %settings=%$settings if(ref($settings) eq "HASH");
  $settings{remove_gaps}=1 if(!defined($settings{remove_gaps}));

  # force these sequences to be just sequences
  my @seq1=@$seq1;
  my @seq2=@$seq2;
  @seq1=map($_->seq,@$seq1) if(ref($$seq1[0])=~/Bio/i);
  @seq2=map($_->seq,@$seq2) if(ref($$seq2[0])=~/Bio/i);

  if($settings{remove_gaps}){
    # find gap positions
    my (@gapPos,%gapPos);
    my $gapString='-';
    my $length=length($seq1[0]);
    for my $seq(@seq1,@seq2){
      for(my $i=0;$i<$length;$i++){
        $gapPos{$i}++ if(substr($seq,$i,1) eq $gapString);
      }
    }

    # remove gap positions
    for(@seq1,@seq2){
      my @seqArr=split(//);
      for(my $i=0;$i<$length;$i++){
        $seqArr[$i]="" if($gapPos{$i});
      }
      $_=join("",@seqArr);
    }
  }

  # give some fake sequence if the group's sequence is missing
  # AKA pseudocount material
  my $randomSeqTemplate=(defined($seq1[0]))?$seq1[0]:$seq2[0];
  if(!$seq2[1]){
    logmsg "Warning: at least two sequences are needed in group 2. Adding missing sequences to group 2.";
    $seq2[0]=join("",LKUtils::Shuffle([split(//,$randomSeqTemplate)])) if(!$seq2[0]);
    $seq2[1]=join("",LKUtils::Shuffle([split(//,$randomSeqTemplate)])) if(!$seq2[1]);
  }
  if(!$seq1[1]){
    logmsg "Warning: at least two sequences are needed in group 1. Adding missing sequences to group 1.";
    $seq1[0]=join("",LKUtils::Shuffle([split(//,$randomSeqTemplate)])) if(!$seq1[0]);
    $seq1[1]=join("",LKUtils::Shuffle([split(//,$randomSeqTemplate)])) if(!$seq1[1]);
  }

  my $numSeq1=scalar(@seq1);
  my $numSeq2=scalar(@seq2);
  $settings{seqs_are_same}=1;
  my $within1=LKUtils::_pDistance(\@seq1,\@seq1,\%settings);
  my $within2=LKUtils::_pDistance(\@seq2,\@seq2,\%settings);
  $settings{seqs_are_same}=0;
  my $between=LKUtils::_pDistance(\@seq1,\@seq2,\%settings);

  # sanity check, but die() shouldn't be invoked here because the sequences should be declared by now
  die "Error: The number of sequences is zero in either group 1 ($numSeq1) or group 2 ($numSeq2) and so Fst cannot be calculated:\n\n".join("\n\n",join("\n",@seq1),join("\n",@seq2))."\n" if(!($numSeq1 && $numSeq2));

  # Normalize by number of comparisons.
  # If there are no comparisons (ie one sequence in a group),
  # then pretend that it was compared with itself.
  # However, there should be enough sequences in the groups at this point and so this is just a sanity check.
  my $within1Comparisons=($numSeq1*($numSeq1-1))/2 || 1;
  my $within2Comparisons=($numSeq2*($numSeq2-1))/2 || 1;
  my $betweenComparisons=$numSeq1*$numSeq2;

  #print "(($within1/$within1Comparisons+$within2/$within2Comparisons)/2)\n";
  my $within=(($within1/$within1Comparisons+$within2/$within2Comparisons)/2);
  $between=$between/$betweenComparisons;

  my $Fst=($between-$within)/$between;

  return $Fst;
}

# p-distance between two groups of sequences
# settings parameter: seqs_are_same: skips distance calculations within the same seq 
sub _pDistance{
  my($seq1,$seq2,$settings)=@_;

  my $alnLength=length($$seq1[0]); # assume flush aln
  my $distance=0;
  for(my $i=0;$i<@$seq1;$i++){
    my $sequence1=lc($$seq1[$i]);
    my $j=0;
    $j=$i+1 if($$settings{seqs_are_same});
    for($j=$j;$j<@$seq2;$j++){
      my $sequence2=lc($$seq2[$j]);
      my $distanceIJ=0;
        $distanceIJ++; #pseudocount
      for(my $k=0;$k<$alnLength;$k++){
        my ($nt1,$nt2)=(substr($sequence1,$k,1),substr($sequence2,$k,1));
        next if(
          $nt1 eq $nt2
          || $nt1 eq '-' || $nt2 eq '-'
        );
        $distanceIJ++;
      }
      $distance+=$distanceIJ;
    }
  }
  return $distance;
}

# reads a shuffled paired end fastq and trims it
# and cleans it
# Author: Lee Katz <lkatz@cdc.gov>
sub qualityTrimFastq($;$){
  my($fastq,$settings)=@_;
  $$settings{numcpus}||=LKUtils::getNumCPUs();
  # trimming
  $$settings{min_quality}||=35;
  $$settings{bases_to_trim}||=10; # max number of bases that will be trimmed on either side
  # cleaning
  $$settings{min_avg_quality}||=30;
  $$settings{min_avg_quality_window}||=70;
  $$settings{min_length}||=62; # twice kmer length sounds good

  my (@t,$t);
  my $Q=Thread::Queue->new;
  for(0..$$settings{numcpus}-1){
    $t[$t++]=threads->create(sub{
      logmsg "Launching thread TID".threads->tid;
      my(@entryOut);
      while(defined(my $entry=$Q->dequeue)){
        my($id1,$seq1,undef,$qual1,$id2,$seq2,undef,$qual2)=split(/\s*\n\s*/,$entry);
        my %read1=(id=>$id1,seq=>$seq1,qual=>$qual1,length=>length($seq1));
        my %read2=(id=>$id2,seq=>$seq2,qual=>$qual2,length=>length($seq2));
        # quality trim: read 1, 3' read
        my($numToTrim3,$numToTrim5,@qual);
        $numToTrim3=$numToTrim5=0;
        @qual=map(ord($_)-33,split(//,$read1{qual}));
        for(my $i=0;$i<$$settings{bases_to_trim};$i++){
          $numToTrim3++ if($qual[$i]<$$settings{min_quality});
          last if($qual[$i]>=$$settings{min_quality});
        }

        $read1{seq}=substr($read1{seq},$numToTrim3,$read1{length}-$numToTrim3-$numToTrim5);
        $read1{qual}=substr($read1{qual},$numToTrim3,$read1{length}-$numToTrim3-$numToTrim5);
        # quality trim: read 1, 5' read
        $numToTrim3=$numToTrim5=0;
        @qual=map(ord($_)-33,split(//,$read2{qual}));
        for(my $i=0;$i<$$settings{bases_to_trim};$i++){
          $numToTrim3++ if($qual[$i]<$$settings{min_quality});
          last if($qual[$i]>=$$settings{min_quality});
        }
        for(my $i=$read2{length}-$$settings{bases_to_trim};$i<@qual;$i++){
          $numToTrim5++ if($qual[$i]<$$settings{min_quality});
          last if($qual[$i]>=$$settings{min_quality});
        }
        $read2{seq}=substr($read2{seq},$numToTrim3,$read2{length}-$numToTrim3-$numToTrim5);
        $read2{qual}=substr($read2{qual},$numToTrim3,$read2{length}-$numToTrim3-$numToTrim5);

        #cleaning stage
        # TODO keep singletons
        next if(length($read1{seq})<$$settings{min_length});
        next if(length($read2{seq})<$$settings{min_length});
        # avg quality
        my $qual_is_good=1;
        QUALITY_CHECK:
        for my $qualStr ($read1{qual},$read2{qual}){
          my $length=length($qualStr)-$$settings{min_avg_quality_window};;
          my @qual=map(ord($_)-33,split(//,$qualStr));
          for(my $i=0;$i<$length;$i++){
            my @subQual=@qual[$i..($$settings{min_avg_quality_window}+$i-1)];
            my $avgQual=sum(@subQual)/$$settings{min_avg_quality_window};
            if($avgQual<$$settings{min_avg_quality}){
              #print "$i: ".join(".",@subQual)."\n$avgQual <=> $$settings{min_avg_quality}\n";die;
              $qual_is_good=0;
              last QUALITY_CHECK;
            }
          }
        }
        next if(!$qual_is_good);

        my $entryOut="$read1{id}\n$read1{seq}\n+\n$read1{qual}\n$read2{id}\n$read2{seq}\n+\n$read2{qual}\n";
        push(@entryOut,$entry);
      }
      return \@entryOut;
    });
  }

  my $entryCount=0;
  open(FQ,'<',$fastq) or die "Could not open $fastq for reading: $!";
  while(my $entry=<FQ>){
    $entry.=<FQ> for(1..7); # 8 lines total for paired end entry
    $Q->enqueue($entry);
    $entryCount++;
    logmsg "Finished loading $entryCount pairs" if($entryCount%100000==0);
    #if($entryCount>99999){
    #  logmsg "DEBUG: only $entryCount loaded";
    #  last;
    #}
  }
  close FQ;
  $Q->enqueue(undef) for(1..$$settings{numcpus});

  # status update
  STATUS_UPDATE:
  while(my $p=$Q->pending){
    logmsg "$p entries left to process";
    # shave off some time if the queue empties while we are sleeping
    for(1..5){
      sleep 2;
      last STATUS_UPDATE if(!$Q->pending);
    }
  }

  my @entry;
  for(@t){
    my $tmp=$_->join;
    push(@entry,@$tmp);
  }

  my $numCleaned=$entryCount-@entry;
  logmsg "$numCleaned entries out of $entryCount removed due to low quality";
  return \@entry;
}


# RANGES
our %minInRange;
our %maxInRange;
sub addNumberForRange{
  my($num,$label)=@_;
  $label="" if(!defined $label);
  LKUtils::max($num,$label);
  LKUtils::min($num,$label);
}
sub max{
  my($num,$label)=@_;
  $label="" if(!defined $label);
  if($num!~/[^\d]/ && (!defined($maxInRange{$label}) || $maxInRange{$label}<$num)){
    $maxInRange{$label}=$num;
  }
  return $maxInRange{$label};
}
sub min{
  my($num,$label)=@_;
  $label="" if(!defined $label);
  if($num!~/[^\d]/ && (!defined($minInRange{$label}) || $minInRange{$label}>$num )){
    $minInRange{$label}=$num;
  }
  return $minInRange{$label};
}
sub range{
  my($label)=@_;
  $label="" if(!defined $label);
  my($min,$max)=(LKUtils::min('x',$label),LKUtils::max('x',$label));
  return ($min,$max) if wantarray;
  return "$min-$max";
}

"Ugh, Perl and its true values";

